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The quantum master equation applied to electronic transport through nanoscopic devices provides 
information not only on the stationary state but also on the dynamics. The dynamics is characterized 
by the eigenvalues of the transition-rate matrix, or generator, of the master equation. We propose to 
use the spectrum of these eigenvalues as a tool for the study of nanoscopic transport. We illustrate 
this idea by analyzing a molecular quantum dot with an electronic orbital coupled to a vibrational 
mode, which shows the Franck-Condon blockade if the coupling is strong. Our approach provides 
complementary information compared to the study of observables in the stationary state. 



I. INTRODUCTION 

Recent progress in nanotechnology allows to fabricate 
transistors with a single molecule forming the active 
rcgionPSEl The transport properties of such devices have 
been studied extensively both experimentally and theo- 
retically. Theoretical approaches have been reviewed by 
Andergassen et aZP and by Zimbovskaya and PedersonP 

The transport properties of single- molecule devices 
are often dominated by strong interactions. Electron- 
electron interaction, which leads to the Coulomb 
blockade}^ and the coupling of electrons to vibrational 
modea^^ have to be taken into account. We assume 
the nuclear motion to be slow on the time scale of 
electronic transitions, which is the case for most, but 
not for all, single-molecule devices studied so farf^ 
The electron-vibron coupling is then due to the change 
of the equilibrium nuclear configuration with the elec- 
tronic state, i.e., the Franck-Condon effect: Even if 
a molecule is initially in a vibrational eigenstate, af- 
ter an electronic transition it will be in a superpo- 
sition of vibrational eigenstates belonging to the new 
electronic configuration. The probability to end up in 
any one of them is determined by Franck-Condon ma- 
trix elements between vibrational eigen states for the old 
and new electronic configurations!^"^ These matrix el- 
ements can be very small if the equilibrium value of 
the relevant normal coordinate changes strongly with 
the electronic state. Consequently, the rates for elec- 
trons tunneling into or out of the molecule and thus the 
current can be strongly suppressed — this is the Franck- 
Condon blockadep3H2i] The dynamics is also unusual in 
this regime: Electrons tunn el in avalanches separated 
by quiet time intervals ! 14 * 15 ! The avalanche-type trans- 
port is self-similar on intermediate timescales and also 

leads to a strong enhancement of the zero- frequency Fano 
factorPHU 

Strong interactions and states far from equilibrium 
make t he description of molecular devices difficult in 
generalP^ We here focus on the case of weak hybridiza- 
tion between the molecule and t he leads. In this case, 
master-equation (ME) approaches^aEMIHUmHH or the 
equivalent real-time diagrammatic approacfpSHSZI can \, e 



employed. In principle, the ME takes into account all 
interactions in the molecule but requires an approximate 
treatment of the coupling to the leads. 

The simplest non-trivial version of the ME treats 
the hybridization to second order (sequential tunnel- 
ing) and neglects off-diagonal components of the re- 
duced densit y matri x in the eigenbasis of the molecular 
Hamiltonian! 13 * 14 -^ The neglect of the off-diagonal com- 
ponents is generally not justified if some of the eigenstates 
are degenerate, since then the choice of eigenstates in the 
degenerate subspaces is arbitrary. In the absence of spin- 
dependent terms in the Hamiltonian, spin degeneracy is 
always present. It is then appropriate to retain all sec- 
ular components of the reduced density operator) 18 ! 28 ! 29 ! 
i.e., all diagonal matrix elements and off-diagonal matrix 
elements between degenerate states. 

Going beyond sequential tunneling, at fourth order one 
obtains cotunneling and pair tunneling as well as a con- 
tribution to the life-time broadening of the molecular 
levelsPSl Koch et aZP^ include cotunneling but consider 
only the diagonal components of the reduced density ma- 
trix. However, integrating out the off-diagonal terms of 
second order generates contributions of fourth order to 
the transition rates between the diagonal components, 
which are taken into account by Leijnse and Wegewijgi 8 -! 
as well as Koller et alF^ 

The dynamics of charge transport through molecules, 
in particular in the Franck-Condon regime, have so far 
mostly been studied by considering the curre nt-noise 
spectrum and the full counting statistics!^ ! 16 * 30 ! For the 
Franck-Condon regime, these studies have been aug- 
mented by real-time Monte Carlo simulationsP^HSD Re- 
cently, Donarini et a/! 3 -^ have considered a harmonically 
varying bias voltage. Assuming spinless electrons and 
employing the Markov and sequential-tunneling approx- 
imations, they find hysteretic behavior of current and 
charge! 3 -^ In the present paper, we propose and illus- 
trate a different approach to the dynamics: We study the 
eigenvalue spectrum of the transition-rate matrix appear- 
ing in the ME. Since these eigenvalues describe the decay 
rates and oscillation frequencies of deviations of the re- 
duced density matrix from the stationary solution, they 
provide a complementary view of the dynamics. More 
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specifically, we will consider the part of the spectrum 
that is small in absolute value, corresponding to states 
that decay or oscillate slowly. We use the well-studied 
molecular transistor with a single vibrational mode as a 
testbed for our idea. 

The remainder of this paper is organized as follows: In 
Sec. [n] the model is briefly introduced and the ME and 
the relevant approximations are discussed. The transi- 
tion-rate matrix is introduced and the physical interpre- 
tation of its eigenvalues is given. Results for the eigen- 
value spectra in various regimes are presented and dis- 
cussed in Sec. Mil Section lPVl sumrnarizes the main results 
and draws some conclusions. 



II. MODEL AND METHOD 

A. Anderson-Holstein Hamiltonian 

In the following we consider a device containing a sin- 
gle molecule coupled to two electrodes. The device is 
described by the Anderson-H olstein Hamiltonian H = 
H leads + H mol + H t CHEESES whe re 



leads / ^\j\<io^v}L<J 



(1) 



refers to the noninteracting leads. For simplicitly, we as- 
sume that both electrodes have identical band structures 
and a constant density of states. The operator cJ, k(T cre- 
ates an electron in lead v = L, R with momentum k and 
spin a. The molecule is described by 



#moi = ^2 £d d ^ dcr + \ ^{hd - 1) 

a 

+ hu v (tfb+ ij +Ahu v (b + tf)n d , (2) 



where d\ creates an electron with spin a and energy e d 
in the molecular orbital, n d = dldf + d\d^ is the corre- 
sponding number operator, and W is the raising operator 
of a harmonic vibration mode. In a break-junction de- 
vice, the on-site energy can be tuned by a gate voltage, 
which is absorbed into e d . Finally, the tunneling between 
the molecule and the leads is described by 



H t = - 



1 



N V 

vvlo 



[tvdl 



H.c. 



(3) 



where N ^> 1 is the number of sites in one lead, which is 
often absorbed into t v W^ We here assume that the tun- 
neling matrix elements t v are independent of momentum 
and spin in the relevant energy range and that the con- 
tacts are symmetric. 



B. Master equation 



The reduced density operator of the molecule is 



Pmoi(i) = Tri cads p(t), 



(4) 



where p(t) is the full density operator of the molecule and 
the leads and the trace is over all many-particle states 
of the leads. The full density operator satisfies the von 
Neumann equation 



dp 

~dt 



(5) 



There are various ways to obtain a ME starting from 
Eq. pSjlP*' Under the condition that the system was in a 
product states with the leads in (separate) thermal equi- 
librium at some initial time to, p(to) = Pmo ifa) < 8> P? ea ds! 
one can derive a ME that is local in timeJEEIBQ] This 
so-called time-convolutionless ME reads 

^" 0l = — jr [-Hmol) Pmol] — *o) Pmol = A p mo l (6) 

where the commutator induces the bare time evolution of 
the decoupled molecule and C(t, t ) is a linear superoper- 
ator describing the coupling to the leads. The right-hand 
side of the ME is linear in p mo i so that we can rewrite it 
as a product involving a superoperator A. 

We write p mo i in the basis of eigenstates \n,q) to the 
eigenvalues E nq of H mo \, where n e {0, t, 4-, t4-} specifics 
the electronic state and q is the quantum number of the 
vibration. The eigenenergies are 



E„ 



u i 
e d n d + — n d (n d - 1) 

+ %uj v (q + ^ - X 2 Tiuj v n 2 d 

(e d - X 2 Huj v ) n d + - (U - 2X 2 Huj v ) n d (n d - 1) 



+ Tiuj v q + 



1 

2/' 



(7) 



where n d (n) = 0, 1,2 is the number of electrons in the 
electronic state n. The matrix elements of the reduced 
density operator in this basis are denoted by 



p nn , 

Hqq> 



(n,q\ p r „oi \n',q'). 
The ME ([6]) written in components reads 
d 



(8) 



— p nn , 
dt' qq 



TP \ n nn 
^n'q'J Pqq 1 



E 

n" ' q"n" 1 ' q"' 



c 



nn ,n n 



n n 

Pq „ ql ,. 



(9) 



Expressed precisely, our goal is to find the eigenvalue 
spectrum of the superoperator A in the ME The 
physical interpretation of the eigenvalues becomes clear 
if one inserts the ansatz 



Pmol(i) 



(10) 
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into Eq. Here, a is a complex number and ( a is an 
operator on the molecular Fock space. This leads to the 
eigenvalue equation 

A( a =a( a . (11) 

The ansatz ( [io| indeed solves the ME if a is an eigen- 
value of A. Then, Re a is the negative of the decay rate 
of the corresponding solution, while Ima is the angu- 
lar frequency of its oscillations. Evidently, a vanishing 
eigenvalue a — corresponds to a stationary state. The 
stationary density operator is thus the right eigenvector 
£o if we impose the normalization condition Tr£o = 1- 

Since the ME ^ has to preserve the trace of p mo i, it 
must satisfy 

J nn 

°=£^f = -£ £ f (12) 

nq nq n" q"n"' q'" 

for all p mo i. Therefore, r] with components Vo!qq' = 
Snn'Sqq 1 is a left eigenvector of A to the eigenvalue zero. 
This proves that at least one stationary state exists. This 
solution is unique if the system is ergodic in the sense 
that every state can be reached from every other state 
by a finite number of transitions.^^ This condition is 
satisfied by our model for non-zero temperature. Thus 
there is exactly one eigenvalue a = 0. 

What is the meaning of the other right eigenvectors 
( a for a ^ 0? These eigenvectors are not well-formed 
density matrices since they have zero trace. This follows 
from the fact that tjq is a left eigenvector to eigenvalue 
zero. Since the left and right eigenvectors to different 
eigenvalues are orthogonal, one has for all right eigenvec- 
tors £ a to non- vanishing eigenvalues 

= TrryJCa= £ «%)* = ( 13 ) 

nqn'q' nq 

Thus one cannot interpret ( a as a density matrix. How- 
ever, linear combinations of the form 

ftnolM = Co + ^c a e a *C Q , ( 14 ) 

with constants c Q chosen such that p mo i is hermitian, 
have unit trace and statisfy the ME ^ . As long as p mo i 
is a positive matrix, it is a permissable time-dependent 
density matrix. Hence, the eigenvectors ( a for a ^ 
describe deviations of possible density matrices from the 
stationary state. Their time dependence is governed by 
the eigenvalue a. One can show that all eigenvalues a^O 
have negative real partsj^S i.e, they describe deviations 
from the stationary state that decay for t — > oo. 

In practice, approximations are needed to obtain the 
superoperator A. We employ the sequential-tunneling 
and secular approximations. While these are standard for 
the study of stationary states, we have to show that they 
are justified for our purpose of obtaining the eigenvalue 
spectrum. We assume weak coupling between molecule 



and leads and expand the right-hand side of Eq. (|6| in 
powers of ^l,_r- For the tunnel Hamiltonian H tl only even 
powers of t L ^ R enter. One can thus write 

j oo 

^ = ]T^) Pmo , (15) 

Here, A>°'p mo i = — (i/K) [Hmo\> Pmoi] is the bare time evo- 
lution from Eqs. (JsJ) and ([9|. Equation Q shows that 
„4(°' is diagonal in the basis {|n, q){n', q'\}. We split the 
reduced density operator into a secular part p s and a 
non-secular part p n and expand both, 

oo 
n=0 

This expansion allows us to write down the ME order by 
order in t\ R . 

First, we discuss the stationary state, for which the 
left-hand side of the ME ^ vanishes. At order zero, we 
obtain 

A^p^+A^p^ =0. (17) 

Since A^ is diagonal and according to Eq. ^ gives zero 
when acting on the secular part p^ , we obtain A^ p n ^ = 
0. This implies that p n °^ — 0, since _4(°) multiplies any 
non-secular component of the reduced density operator 
by a non- vanishing factor. At second order we thus find 

A (?) p (S» + A (2) p (0) + A Q) p W + A (0) p (2) 

= AWpW+AW p ™ =0. (18) 

The leading secular components are thus of order zero, 
while the non-secular components are at least of order 
two. From the ME up to second order in t^ ^, one cannot 
obtain the second-order corrections to the secular part, 

p^s \ this requires to go to fourth order P^HHxhe leading- 
order stationary density operator is thus the solution of 

A^p^ = (19) 

with the non-secular components vanishing, pn^ = 0. 

Next, we turn to the dynamics. In the ME the 
time derivative of all non-secular components of /3 mo i 
picks up imaginary factors —{i/Ti) (E nq — E n i q i). They 
are large compared to the typical scale introduced by 
the coupling since we assume the coupling to be weak 
and exclude accidental near-degeneracies. To put it dif- 
ferently, the zero-order superoperator A^ ' has eigen- 
values a q ° q ]' nn = —(i/K) (E nq — E n > q >) to eigenstates 
\n, q)(nf, q'\. These eigenvalues are zero for secular com- 
ponents and large and purely imaginary for non-secular 
components. As long as the perturbative expansion in 
tL,R is justified, the full superoperator A will have eigen- 
values close to the zero-order ones, thus some will be 
small in absolute value, while the rest will have a large 
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imaginary part. We are interested in the part of the spec- 
trum with small absolute value. This is the part lacking 
large imaginary parts from order zero. Consequently, to 
find the small eigenvalues, one has to consider the secular 
sector. To leading order, the eigenvalues arc then given 
by the second-order superoperator A^ . Consequently, 
we have to solve the ME 



dp s 
dt 



= A {2 



(20) 



for the secular density operator p s . T he derivation of 
A&> is standard and we omit the details.' 14 ' 16 ' 22 ' 27 ' 28 ' 46 ' 47 ' 
Assuming that at some initial time t the molecule and 
the leads are in a product state with the leads in separate 
equilibrium, p(t ) — p m oi(to)®Pi ca ds> anc ^ taking the limit 
to — > — oo, one obtains to second order in ti, u, 

dpmol i rjr 

— — r L"mol)PmolWj 



1 

h 2 j,, 



drTt 



loads 



[Ht, [e 



-»(ffload B + -Hmol)T/R 



x H t e « H >~*+ B «»*W\ p mol (t) ® p? oads ]] • (21) 
We expand the nested commutators and take the trace 



over the lead degrees of freedom using 



Tri ca ds Pk Ja ds c tkcr C ' / 'k'CT' — <W'<^kk'<J<rcr'/(£r/k)! 



(22) 



Trieads Pl ea ,ds c vk<T c l'k><T< = <W<W<W[ 1_ /(^k)], (23) 

where £„k = £k — p v , p v is the chemical potential in 
lead v, and /(£) is the Fermi function. The chemical 
potentials in the left and right leads satisfy pr — Pl — 
eV. We assume that the device is symmetric, i.e., ph — 
~Pr — —eV/2. In the basis of eigenstates \n, q), the ME 
has the form 



d o nn 

" aa ' C 1 "\ n nn ' — \ , nnn" n"n' 

— ^ \J^nq J^n'q' ) P qq ' / y "•qq" Pq"q' 

n" q" 



dt 



E/ V "'..'•/ 

n" q" 



( r>n'n"\* nn" , 
[Rq'q" ) Pqq" + 



E 

n"n nr q" q'" 



r>nn" ,n"'n' n"n"' 
qq" \q"'q' Pq"q"' ' 



(24) 



Since we are only interested in the secular sector, we 
can drop the first term on the right-hand side, which 
corresponds to A 1 -- * 1 . The expressions for the rates R also 
simplify in this sector. The rates appropriate for secular 
p mo \ read 



nnn 
n qq' 



R 



nn ,n n 
qq" ,q"'q> 



E 



n" q" 



eV 



eV 



7 2 J \**'n"q" ^nq \ 2 

^ nn" n" n' r qq r q" q' ~ / , nn" ,±J n" n' r qq" r q" q 
(7 a 

(^E n q ~ E n iiq" — ^ + / \E nq — E n 'i q ii + — J 

En" „D t<T , F »F f + Y^Z? tcr D° ' ,,, ,F^ F m 
1J nn"^n'"n' ± qq 1 q"'q' ^ / , u nn""n'"n" qq" 1 q q 



(25) 



(26) 



where the rate T = 2TrNL y n\tL,R\ 2 /ft describes the cou- 
pling to the electrodes with densities of states Nl,r, 
which are assumed to be constant, and DZ„> = (n\d a \n'), 



D 



to- - 



A4\ 



are matrix elements of the electronic 



operators. The Franck-Condon matrix elements F qq / = 
read explicitly^M 1 ^ 



(e | e -A(6t-6)| 9 r 



F i — 
r qq ~ 



\i>-9< e - x / 2 i*>-*<(A 2 ) 



1 



for q> q' , 
for q < q' , 



(27) 



and F 



,„, - {Fq'qY = F q'q- Her e, q< = mm(q,q{), 
g> = max(g, q'), and L^(x) are generalized Laguerre 
polynomials. 

and 



D ° n , always appear in combinations corresponding to the 
creation and annihilation of electrons of the same spin 
(j. This leads to the vanishing of the rates for certain 
combinations of electronic states. Furthermore, the only 
off-diagonal secular components of p mo \ are pj^ and p^J 
for all q. Thus, the rates relevant for the secular sector 
simplify to 



In Eqs. (25) and (26), the matrix elements D° 
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r>nn 

qq 



n qq",q"q 



n" q" 



eV 



f ( E n ,> q >> - E nq — J + / 1 E n n q „ - E nq + t) 



a (7 

eV 



Ei^»"i 2 i^"i 2 +Ei^"«i 2 i F «" 



eV 



.1 [ E nq — E n n q n — J + / ( E nq - E n "q" + — 



(28) 



(29) 



r 



Note that principal-value integrals, which plague the ME 
for the full reduced density operator, cancel in the secular 
sector. Also, the rates in Eqs. ( 28 ) and ( 29 1 are real. The 



diagonal components of the ME ( 24 ) then simplify to 



dp. 



,nn 

qq 



dt 



~ L ^qq Pqq 



Ep»»" ,n"n n"n" l nn\ 

H qq",q"q Pq"q" > W) 



n" q" 



while for the off-diagonal secular components we obtain 



The current requires a different approach. In the follow- 
ing, we consider a charge current flowing from left to right 
as positive. The operator I v for the current between the 
lead v — L,R and the molecule is 



d A 
ve — N v 

dt 



ve - [H, N v ] 
n 



i^[HZ,ft v ), (32) 



dp, 



dt 



nn 



r>n n 
^qq 



)P. 



a. n 

qq 



(31) 



The equations for the diagonal and the off-diagonal com- 
ponents thus decouple and the off-diagonal components 
exhibit simple exponential decays. 

Averages of local observables such as the electron num- 
ber in the molecule can be obtained directly from p qq , . 



where the numerical value of v is +1 (—1) for the left 
(right) lead, — e is the charge of the electron, N v is the 
number operator of electrons in lead is, and fl" is the 
part of the tunneling Hamiltonian H t involving lead v. 
The current is then l v = Tr I v p, where Tr is the trace 
over the full Fock space. Under the same assumptions as 
used for the rates above one obtains 



i v = vev y: e 



nqn' q' n" q 
— f\ E n " q " — E n 



f[E n , 



E nq — V ■ 



eV 



n'n" n"n M q'q" c . 



q"q 



eV 



1J n'n" 1J n"n F q q r q"q 



p nn , 
Pqq' 



(33) 



where the sum over n, q, n', q' only runs over secular 
components. 

For the numerical calculations, we cut off the ladder 
of harmonic-oscillator states so that < q < q max . Then 
the dimension of the molecular Fock space is 4(g max + 1). 
The secular part of p mo \ has 6(g max + 1) components — 
the 4(g max + l) diagonal ones and 2(g max + f) off-diagonal 
ones of the form p^jj and p^. In the secular sector, the 
density operator can thus be represented by a 6(g max + I)- 
component vector and the superoperator by a real 
6(<Zmax + 1) x 6(g max + 1) transition-rate matrix. 

It is a common problem of the ME approach that the 
matrix can be ill-conditioned, i.e., the ratio of the 



largest to the smallest eigenvalue can be very large. This 
is here due to the Fermi functions and Franck-Condon 



matrix elements in Eqs. (28) and (29), which span many 



orders of magnitude. For this reason, black-box diagonal- 
ization routines often fail to distinguish the true station- 
ary state from eigenvectors to very small eigenvalues. To 
overcome this difficulty, we use Mathematical to solve 
the eigenvalue problem with high precision. We adapt 
the number of digits in the calculation such that it is 
larger than the L°° condition number of the matrix by 
at least 12, which should give results for the eigenvalues 
and eigenvectors with on the order of 12 significant digits. 

Before analyzing eigenvalue spectra in the next section, 



(> 



we comment on their dependence on the cutoff <7 max - We 
find that adding highiy excited vibrational states only 
adds eigenvalues to the middle and upper part of the 
spectrum. If g max is not too small, the spectrum at small 
magnitudes does not change significantly with (? max . This 
is plausible since adding highly excited vibrational states 
should not introduce additional slow relaxation channels. 
We use (? m ax = 30, unless noted otherwise. 

It would be possible extend the analysis to higher or- 
ders in the tunneling amplitudes th,R- This would re- 



quire to solve the eigenvalue equation ( 11 ) perturbatively. 



Compared to time-independent perturbation theory for 
Hamiltonian systems, this is more complicated since the 
superoperator A is not hermitianPsl In this work, we ob- 
tain the eigenvalues a up to second order. The next non- 
vanishing contribution is of fourth order. At this order, 
the restriction to the secular sector is not possible, we re- 
quire A^ for all secular and non-secular states, and we 
also need A^ within the secular sector. The s tation ary 
solution at this order has been studied beforeP32U We 
leave the spectral analysis for a future work. 



III. RESULTS 

In the present section, we present results for the eigen- 
value spectrum in the regimes of the transmitting quan- 
tum dot, the Coulomb blockade, and the Franck-Condon 
blockade. It is shown that the spectrum differs qual- 
itatively between these cases and is characteristic for 
each. We compare the information that can be gained 
more conventionally from observables such as current and 
charge in the stationary state. 




I L_ 
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FIG. 2: (a) Real and (b) imaginary parts of the eigenvalues 
of the transition-rate matrix as functions of the applied bias 
voltage eV. The parameters are the same as in Fig. [T] 




-15 -10 -5 5 10 



bias voltage sV/to, 

FIG. 1: (Color online) (a) Current, (b) electronic occupation 
number, and (c) excitation of the harmonic oscillator as func- 
tions of bias voltage eV for relatively small electron-vibron 
coupling A = 1, on-site energy e<i = 1, Hubbard interaction 
U = 6, and thermal energy fcsT = 0.05. All energies are 
given in units of the vibron energy fiu v — 1. 



A. Transmitting regime 

First, we consider the transmitting molecular device. 
We tune the effective on-site energy — \ 2 hu> v in Eq. Q 
to zero so that resonant tunneling is possible at vanish- 
ing bias voltage. Figure [T] shows the current II, the elec- 
tronic occupation number (n^), and the vibron excitation 
(q) as functions of the bias voltage in the stationary state 
for relatively small electron-vibron coupli ng A = 1 . The 
current increases with characteristic stepsP^EHEZl There 
is a step at zero bias since the device is on resonance. 
The steps at non-zero bias result from inelast ic tunn el- 
ing under excitation of 1,2,... vibron quanta,^ 3 * 14 * 22 ! in 
agreement with the observed average excitation (q) . The 
average electronic occupation changes only weakly and 
mostly above an energy scale on the order of U. This 
weak dependence is due to a small admixture of doubly 
occupied states at higher bias voltages. 

The real and imaginary parts of the eigenvalues a of 
the transition-rate matrix are shown in Fig. [2] for the 
same parameters. The eigenvalue zero is always present, 
as it must be. The real and imaginary parts of the other 
eigenvalues reflect the step positions from Fig. [T] except 
for the zero-bias step. Note that at any voltage, most 
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eigenvalues have vanishing imaginary parts. The eigen- 
values with non- vanishing imaginary part form complex- 
conjugate pairs, since the transition-rate matrix is real. 
The non-vanishing imaginary parts are typically small 
compared to the real parts. This means that the decay 
time of the corresponding deviations from the stationary 
state is much shorter than their oscillation period. Simi- 
lar behavior is found for randomly distributed transition 
rates, where it is essentially a consequence of different 
scaling of the real and imaginar y p arts with the dimen- 
sion of the molecular Fock spaced 

At zero bias, all eigenvalues are real. This has to be 
the case since for V = the molecule is coupled to 
an equilibrium bath and all transition rates satisfy de- 
tailed balancePSES This is easily confirmed by check- 



ing that the non- vanishing rates in Eq. ( 30 1 satisfy 

r,n n,nn / -r-,nn ,n n (E„„—E 



,)/kBT at V = 0. Then 



i <i,m ' n .9 9 
for the diagonal components, the transition-rate matrix 

A can be written in the form 



■Aij — 



for i ^ j, 



(34) 



where j3 = 1/fcgT is the inverse temperature and R®j = 
R ^. Introducing the diagonal superoperator O with com- 
ponents Oij = 8ij e /3Ei / 2 , one easily sees that OAO~ x is 
real and symmetric and therefore has only real eigen- 
values. Since this is a si milarity transformation, A has 
the same real eigenvalues-EMED The 

previous argument 

only applies to the diagonal components. However, the 
off-diagonal secular components show a simple exponen- 
tial decay anyway, as expressed by Eq. ( [31] ). 

Concerning our goal of characterizing different regimes 
in terms of their eigenvalue spectra, the crucial observa- 
tion is that the spectrum shows a clear gap in the real 
part. Thus there are no slow modes — all deviations from 
the stationary state decay with rates that are on the or- 
der of the characteristic rate T. 

It is interesting to analyze the character of the station- 
ary state and of the deviations that decay most slowly. 
At V = 0, the stationary (equilibrium) state is a mix- 
ture of the microstates \n,q) = |0,0), |t, 0), ||,0) with 
equal probabilities, except for exponentially small ther- 
mal occupations of higher-g and doubly occupied states, 
see Figs.[l|b) and[TJc). 

As discussed, the eigenvectors to non-vanishing eigen- 
values represent deviations from the stationary state. 
The components Ca U qq with the largest magnitudes char- 
acterize the microstates that have the largest weight in 
a given deviation. At V = 0, we find that the eigenvalue 
with the smallest non-vanishing magnitude is actually 
threefold degenerate — it corresponds to three linearly in- 
dependent deviations. The corresponding subspace is 
spanned by the hermitian matrices |f, 0)(t, 0| — 1|, 0) (|, 0|, 
|t,0)U,0| + ||,0)(t,0| and -i|t,0)U,0| + i ||, 0)(t, 0|. 
These matrices can be written as 



where a s are the Pauli matrices and |0)(0| is the projec- 
tion operator onto the q — vibron state. The expres- 



sion ( 35 ) shows that the slowest deviations represent spin 
polarizations in the x, y, and z direction in the singly oc- 
cupied sector. Evidently, spin polarizations decay most 
slowly. We will return to this point below. 




FIG. 3: Secular components Ca^<j °f ( a ) the stationary density 
matrix and (b) the slowest deviation for bias voltage eV = 3 
in units of the vibron energy Hlu v — 1. The other parameters 
are the same as in Figs. [I] and [2] 
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FIG. 4: Secular components Ca™ qq of (a) the stationary density 
matrix and (b) the slowest deviation for bias voltage eV = 11 
in units of the vibron energy hu) v — 1. The other parameters 
are the same as in Figs. [I] and [2] 

At the first step in Fig.[T] where eV/2 ps Huj v , there is a 
crossing in the spectrum in Fig. [2](a) and we thus expect 
the deviation with the smallest decay rate to change in 
character. Figure [3] shows the secular components of the 
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stationary density matrix and of the slowest deviation at 
eV — 3fruj v . Compared to V = 0, the stationary state 
obtains finite probabilities for low-lying vibron excitions 
in the sectors of electronic occupation numbers and 1. 
The deviation with the smallest decay rate is now non- 
degenerate and the large components C"" g change sign 
when the occupation changes between and 1 and also 
when q is increased by unity. How can we understand 
this? At eV — 3 Tiuj v , q can increase at most by unity in 
a sequential-tunneling event. Sequential tunneling also 
changes the occupation by ±1. The slowest deviation 
is dominated by an imbalance between the probabilities 
of the microstates |0, 0) , |t, 1), 11,1), |0, 2) on the one 
hand and of |t,0), ||,0), |0,1), |f,2) s \±,2) on the other. 
This imbalance relaxes slowly because endothermal tran- 
sitions between any microstate from one class and any 
microstate from the other are thermally suppressed. 

At the third step at eV « 6 hco v , there is another cross- 
ing in Fig. [5Ja) . Figure [4] shows the secular components 
of the stationary density matrix and of the slowest devi- 
ation at eV = 11 huj v . The stationary state now contains 
highly excited vibrons. Also, the probabilities of dou- 
bly occupied states are comparable to those of empty 
and singly occupied states. The slowest deviation is non- 
degenerate and mainly involves a transfer of weight be- 
tween weakly excited vibron states with q < 5 and highly 
excited states with q > 5. The significance of the num- 
ber of 5 becomes clear by inspecting the spectra in Fig. 
[2] At eV = 11 huj v , q can increase by at most 5 in a 
sequential-tunneling event, whereas any decrease is pos- 
sible. The deviation sketched in Fig. |4|b) is mainly an 
imbalance between vibron states that differ in q by more 
than 5. Such a deviation is slow to relax by endothermal 
sequential tunneling since the relaxation requires more 
than one transition. We have checked this interpretation 
by following the slowest deviation to higher voltages. It 
retains its character but the zero of C^ qq shifts to higher 
q (not shown). This is expected since for higher voltages 
larger changes of q in a single transition become possible. 
Figure [2]ja) shows that the slowest modes become slower 
with increasing voltage, due to the decrease of Franck- 
Condon matrix elements F qq i for larger \q — q'\. 



B. Coulomb blockade 

If the molecular energy level td is detuned from res- 
onance, there is a non-zero excitation energy between 
states with different occupation numbers and the current 
is suppressed at small bias voltages. First, we consider 
the case that the molecular energy level td lies below the 
Fermi energy of the leads at zero bias but the addition 
energy e d + U for a second electron lies above. In this 
Coulomb-blockade regime, the current is suppressed by 
the Coulomb repulsion U. For this regime, Fig. [5] shows 
the current II, the electronic occupation (rid), and the vi- 
bron excitation (q) as functions of bias voltage in the sta- 
tionary state for relatively small electron- vibron coupling 




-5 5 

bias voltage eVWco 



15 



FIG. 5: (Color online) (a) Current, (b) electronic occupation, 
and (c) vibron excitation as functions of bias voltage eV for 
relatively small electron- vibron coupling A = 1, on-site en- 
ergy ed = 0, Hubbard interaction U = 6, and thermal energy 
ksT — 0.05. All energies are given in units of the vibron 
energy hoj v = 1. 



A = 1. Compared to the transmitting regime, we observe 
Coulomb blockade for small bias voltages (\eV\ < 2hu v ), 
where all three observables are approximately constant 
and the stationary state is an equal mixture of the degen- 
erate singly occupied ground states |f, 0) and ||, 0) with 
exponentially small corrections. When the bias voltage 
reaches a certain threshold, electrons can tunnel out of 
the molecule so that the average occupation number de- 
creases, see Fig. [5|b), and the current sets in, see Fig. 
[5ja). For higher voltages, also doubly occupied states 
occur with significant probability and the average occu- 
pation number increases again. Similarly to the transmit- 
ting regime, the current increases in steps due to inelastic 
tunneling under excitation of vibrons, as seen from the 
increase in (q) in Fig. [5jc). 

Figure [6] shows the real and imaginary parts of the 
eigenvalues for the same parameters used in Fig. [5] 
Clearly, as the system enters the Coulomb blockade, a 
non-zero real eigenvalue becomes very small. This is a 
threefold degenerate eigenvalue corresponding to devia- 
tions of the form (35). Thus in the Coulomb blockade, 



the spin polarization in the singly occupied sector decays 
very slowly. This is easy to understand: An electron has 
to tunnel out of the molecule and another electron with 
opposite spin has to tunnel in (or vice versa) to relax the 
spin. But the first tunneling process is thermally sup- 
pressed by the exponentially small tail of the Fermi func- 
tion. If we were to include higher orders in t^ R in the 
calculation, eigenvalues with exponentially small leading- 
order contribution would generically obtain a contribu- 
tion of order |£l,.r| 4 that is not exponentially suppressed 
but is still small as long as the pertubative expansion in 
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bias voltage eVf~h(a 



FIG. 6: (a) Real and (b) imaginary parts of the eigenvalues 
of the transition-rate matrix as functions of the applied bias 
voltage eV. The parameters are the same as in Fig. [5] 



tL,R is justified. We have shown above that the same spin 
deviations still decay slowly, although not with exponen- 
tially suppressed rate, in the transmitting regime. Note 
that the inclusion of the off-diagonal secular components 
of Pmoi is necessary to obtain the correct spin symmetry 
and degeneracy of these slow modes. 

Next, we turn to the two regimes where both ed and 
ed + U lie either above or below the Fermi energy of the 
leads. Then, the molecular orbital in the stationary state 
is either predominantly empty or doubly occupied, re- 
spectively. The two regimes are related to each other 
by a particlc-hole transformation so that the transport 
properties are very similar. In these regimes it is the 
single-particle energy rather than the Coulomb interac- 
tion that suppresses sequential tunneling. We neverthe- 
less continue to use the term "Coulomb blockade" . For a 
predominantly empty or doubly occupied molecular or- 
bital, the molecular spin is essentially zero and its relax- 
ation should not be important for the dynamics, like it 
was in the previous case. We plot the stationary current, 
electronic occupation, and vibron excitation as functions 
of the bias voltage for = 4.5 fnu v in Fig. [7j There is 
now a broad regime at low bias voltage where the molec- 
ular orbital is essentially empty. Singly occupied states 
become available above the Coulomb-blockade threshold 
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FIG. 7: (Color online) (a) Current, (b) electronic occupation, 
and (c) vibron excitation as functions of bias voltage eV for 
relatively small electron- vibron coupling A = 1, on-site energy 
ed — 4.5, Hubbard interaction U — 6, and thermal energy 
ksT — 0.05. All energies are given in units of the vibron 
energy hoj v = 1. 



so that a sequential-tunneling current sets in. Vibrons 
start to be excited at the same point since an electron 
tunneling out of the molecule has sufficient excess energy 
to excite the vibration. The corresponding eigenvalue 
spectra are plotted in Fig. [8j It is striking that in this 
case no eigenvalue becomes small right at the thresh- 
old at eV rs ±7 Tilo v — the gap in the spectrum persists 
into the Coulomb-blockade regime. An eigenvalue ap- 
proaches zero, i.e., the gap closes, only at a voltage of 
eV ~ ±5huj v . For smaller voltages, even deeper in the 
Coulomb-blockade regime, there are additional transi- 
tions where further eigenvalues become small. Note that 
the stationary observables in Fig. [7] are all exponentially 
suppressed here. 

The stationary state is of course dominated by |0, 0) 
throughout the blockade regime. We now analyze the 
deviations that become slow as the voltage is lowered. 
At eV sa ±5Huj v , a non-degenerate mode becomes slow 
that mainly involves transfer of weight between |0, 0) 
and excited vibrational, and to a lesser extend electronic, 
states. This slow mode is sketched in Fig. [^a). Below 
the voltage eV ps ±5Tluj v , the excited-state-to-excited- 
state transitions from |0, g + to \f, q) and \i,q) become 
suppressed. In particular, the only rapid decay channel 
of the state |0, 1) (to |t, 0) and ||,0) and then to |0, 0) ) 
is suppressed. Therefore, the slowest deviation mostly 
involves transfer of weight between |0, 1) and |0,0). 

Next, at eV f=a ±3tiui v , another non-degenerate mode 
becomes slow. It is sketched in Fig. |9jb). This mode 
involves a transfer of weight between the two lowest vi- 
brational states |0,0), |0, 1) on the one hand and mainly 
the next state |0, 2) on the other. It becomes slow because 
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FIG. 8: (a) Real and (b) imaginary parts of the eigenvalues 
of the transition-rate matrix as functions of the applied bias 
voltage eV. The parameters are the same as in Fig. [7] 



below this voltage also the decay of |0, 2) is suppressed. 
Analogously, at eV ~ huj v , a further mode sketched in 
Fig. [9jc) becomes slow due to the suppression of the de- 
cay of |0, 3). If we would increase the on-site energy fur- 
ther by means of a gate voltage, we expect more and 
more slow modes to appear. 



C. Franck-Condon blockade 

In this subsection, we turn to the signatures of Franck- 
Condon blockade in the spectra. Like for the trans- 
mitting regime, we tune the effective on-site energy 
e c i — X 2 huj v to zero. Then resonant tunneling is possible 
at V = and any suppression is due to Franck-Condon 
blockade. Figure [T0| shows the stationary current, elec- 
tronic occupation, and vibron excitation as functions of 
the bias voltage for intermediate electron-vibron cou- 
pling A = 2. We choose a larger Hubbard interac- 
tion U = 12 Tilo v since for the previously used value of 
U = 6 fito v , the effective interaction in Eq. ([7| would be- 
come attractive. The main effect of the stronger electron- 
vibron coupling is the suppression of the zero-bias current 
~7aY 
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FIG. 9: Secular components CS^q of deviations becoming slow 
within the Coulomb-blockade regime in Figs. [7] and [8] Panels 
(a), (b), (c) show the modes becoming slow at eV — 5hui v , 
3huj v , huj v , respectively. The other parameters are the same 
as m Figs. [7] and [8] 
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FIG. 10: (Color online) (a) Current, (b) electronic occupa- 
tion, and (c) vibron excitation as functions of bias voltage 
eV for intermediate electron-vibron coupling A = 2, on-site 
energy — 4, Hubbard interaction U = 12, and thermal en- 
ergy fcflT = 0.05. All energies are given in units of the vibron 
energy Tiuj v — 1. 



The corresponding eigenvalue spectra are plotted in 



Fig. 11 We find a smaller gap at low bias voltage, com- 
pared to the case of A = 1 shown in Fig. |2| At V = 0, 
the smallest eigenvalue is threefold degenerate and cor- 
responds to spin imbalances of the form ([35]). The slow- 
est deviations are thus the same as for the transmitting 
regime, but their decay rate has become even smaller. At 
first glance, it might be surprising that the enhancement 
of electron-vibron coupling leads to suppressed spin re- 
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FIG. 11: (a) Real and (b) imaginary parts of the eigenvalues 
of the transition-rate matrix as functions of the applied bias 
voltage eV. The parameters are the same as in Fig. |10| In 
particular, the electron-vibron coupling is A = 2. 



laxation. The reason is that in order for the spin to relax, 
electrons have to tunnel in and out of the molecule. At 
low voltage, the only available transitions are between 
|0,0) on the one hand and |f,0) and ||, 0) on the other. 
But these transitions are now suppressed by the small 
Franck-Condon matrix element -Foo = e ~ X ^ ■ The next 
eigenvalue, which is comparable in magnitude, is not de- 
generate and corresponds to an imbalance between the 
empty and singly occupied states. It becomes slow for 
the same reason. 

Finally, we turn to the case of even stronger electron- 
vibron coupling A. The effective on-site energy e^ — \ 2 huj v 
is again tuned to zero. 



Figure 12 shows the stationary 
current, electronic occupation, and vibron excitation as 
functions of the bias voltage for strong electron-vibron 
coupling A — 4 and U — 40 hoj v . Due to the large value of 
U, a larger cutoff q max = 50 is chosen here. Note the cur- 
rent scale in Fig. 12 [&): The current is strongly reduced 



in magnitude for all voltages, in particular for small ones, 
by the Franck-Condon blockadeP^KOD i n this regime, the 
voltage dependence of the occupation number and of the 
vibron excitation are also suppressed. The corresponding 
eigenvalue spectra are plotted in Fig. [13] The typical real 
and imaginary parts have become smaller and the gap is 
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FIG. 12: (Color online) (a) Current, (b) electronic occupa- 
tion, and (c) vibron excitation as functions of bias voltage 
eV for strong electron-vibron coupling A = 4, on-site energy 
ed = 16, Hubbard interaction U = 40, and thermal energy 
ksT — 0.05. All energies are given in units of the vibron 
energy hu) v = 1. 



completely filled in at all bias voltages shown here. Thus 
there are slow modes in the whole voltage range. At least 
at small voltages, the character of the slowest modes is 
the same as for A = 2, though: The most long-lived devi- 
ations are spin and charge imbalances, the decay of which 
is suppressed by small Franck-Condon matrix elements. 
The inset in Fig. 



12 shows details on the small real 



parts on a logarithmic scale. The small real parts roughly 
follow a log-uniform distribution for a certain range of 
rates. Within this range, the probability density func- 
tion is approximately P(|Rea|) ~ l/|Rea|, i.e., it is 
scale-invariant. The uniform distribution of ln|Rea| is 
caused by the approximately exponential dependence of 
the Franck-Condon matrix elements F qq i on q and q' 
for q,q' <C A 2 . The approximate scale-invariance of 
the distribution of small rates implies that the dynam- 
ics of the system within a certain time window is also 
scale-invariant. This is consistent with the approximate 
self-similarity of the tim e-dep endent transport found by 
Monte Carlo simulations P 11 ^ 1 



IV. SUMMARY AND CONCLUSIONS 

In the present paper, we have studied the eigenvalue 
spectrum of the transition-rate matrix in the ME for a 
molecular quantum dot coupled to metallic leads. The re- 
laxational and oscillatory dynamics of deviations of the 
system from the stationary state are characterized by 
the real and imaginary parts of these eigenvalues, respec- 
tively. We have mainly considered the small eigenvalues, 
which describe the slow dynamics. Conceptually, this is 
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FIG. 13: (a) Real and (b) imaginary parts of the eigenval- 
ues of the transition-rate matrix as functions of the applied 
bias voltage eV . Inset: Absolute value of the smallest non- 
vanishing eigenvalues on a logarithmic scale. The parameters 
are the same as in Fig. |12| In particular, the electron-vibron 
coupling is A = 4. 

similar to analyzing the spectrum of low-lying eigenen- 
ergies of a Hamiltonian system. We have applied this 
idea to a molecular transistor with an electronic orbital 
coupled to a vibrational mode. 

The spectra differ qualitatively between a transmitting 
device, a molecule in the Coulomb-blockade regime, and 
a molecule in the Franck-Condon-blockade regime. We 
demonstrate that the character of deviations from the 
stationary state can be analyzed by considering the large 
components of the corresponding eigenvectors. Some of 
the deviations with the smallest decay rates represent 
non-zero spin polarizations of the molecule. They occur 



in groups of three degenerate modes corresponding to 
polarizations in the x, y, and z direction. In order to ob- 
tain these modes, all secular components of the reduced 
density matrix have to be included. 

In the transmitting regime, the spectrum has a gap 
for any bias voltage, i.e., there are no slowly decay- 
ing deviations on the scale of the sequential-tunneling 
rate T. In the Coulomb-blockade regime with predom- 
inantly single occupation of the molecular orbital, the 
gap in the spectrum closes since relaxation of the elec- 
tronic spin becomes slow. If instead the molecular or- 
bital is predominantly empty or doubly occupied, there 
is no finite spin polarization and thus these slow modes 
do not exist. In these cases, the gap persists into the 
Coulomb-blockade regime. However, deep within these 
regimes the gap closes and more and more modes be- 
come slow at consecutive steps. These modes become 
slow since excited-state-to-excited-state transitions are 
thermally suppressed. The dynamics here contains ad- 
ditional information not accessible by observables in the 
stationary state, which show an exponentially suppressed 
voltage dependence. 

For stronger electron-vibron couping we find that the 
gap becomes small even if resonant tunneling is ener- 
getically possible, since certain transition rates are sup- 
pressed by small Franck-Condon matrix elements. In the 
strong Franck-Condon-blockade regime, the gap closes 
over a broad range of bias voltages since many devia- 
tions now decay slowly. We also find an approximately 
scale-invariant distribution of the slowest rates, consis- 
tent with t he pre viously observed self-similar dynamics 
in real timeP^^ 

In the present paper, we have concentrated on the sta- 
tionary and long-lived states. The spectra obtained in 
this paper show additional structure that we have not 
discussed, suggesting that much more information can 
be extracted from the spectra and the eigenmodes. 
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